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Abstract In this paper, we propose a new method to forecast the drift of objects in near coastal 
ocean on a period of several weeks. The proposed approach consists in estimating the probability 
of events linked to the drift using Monte Carlo simulations. It couples an averaging method which 
permits to decrease the computational cost and a statistical method in order to take into account 
the variability of meteorological loading factors. 
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1 Introduction 

Drift of things in the ocean is potentially dangerous for human activities and marine ecosystems. 
For instance, drifting containers may cause serious accidents in the event of collision with ships and 
oil spills may have very negative impacts especially in coastal areas. In this paper, we investigate a 
new method to forecast the drift of an object in near coastal ocean on a period of several weeks. 

The motion of a drifting object on the sea surface is the net result of a number of forces acting 
upon it (water currents due to tide wave, atmospheric wind, wave motion, wave induced currents, 
gravitational force and buoyancy force). It is possible to estimate the drift trajectory given infor- 
mation on the local wind, the surface current, and the shape and the buoyancy of the object. For 
instance, in order to estimate the position of lost containers, the safety-and-rescue services generally 
use short-term meteorological forecasts as forcing of an hydrodynamic model of drift (see Daniel et 
al. [5]). It is usual in such problems to consider several possible buoyancy and drift properties for 
the object since these features are not known precisely in most cases. An uncertainty on the initial 
conditions (position and time) may also be taken into account (see Hackett et al. |13j). 

Since we are interested in longer periods of time in the present study, we cannot directly use 
meteorological forecasts to estimate the object trajectory. Then the drift forecast has to be led out 
in terms of probability. This will be done using a Monte Carlo method, which makes it possible 
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to estimate the probability of some scenarios linked to the object's trajectory, like the probability 
of being in a given point at a given time or the probability of running aground in given areas, 
for example. It consists in computing the object trajectories corresponding to a large number of 
meteorological time series representative of the climatology on the considered area. Because of the 
variability of the meteorological conditions, it is necessary to compute a large number of trajectories 
in order to get reliable estimates of the quantities of interest. Hence we are faced with two problems. 

At first, as the existing data sets describe the meteorological conditions only on the few last 
decades, it is necessary to be able to simulate new realistic meteorological time-series. For this, we 
have used stochastic models. More precisely, it is assumed that these time series can be decomposed 
as the sum of two components, the first one that describes the meteorological conditions at a synoptic 
scale and the second one that represents fluctuations at the smaller scale. Then, the synoptic 
component is simulated with a non-parametric resampling algorithm proposed by Monbet, Ailliot 
and Prevosto [T^, and the short-term component with an Autoregressive model. 

Secondly, we have to compute the corresponding object trajectories, which oscillate with tide, 
in a reasonable computational time. Indeed, the trajectory of a drifting object submitted to tide 
wave currents and currents generated by other meteorological factors, is essentially composed by a 
trend and oscillations due to the tide, which have a small period with respect to the time period of 
interest (several weeks). Due to these high frequency oscillations, the numerical integration of the 
considered system is generally time consuming. So that it is convenient to split the estimation of the 
displacement of the object in two steps: first, we calculate the trend which is of low computational 
cost thanks to the lack of oscillations and then we reconstruct the oscillation around the trend. 
For this, we apply the Averaging Method developed by Frenod [7] in order to identify the averaged 
fields governing the trend and the oscillating operators allowing the reconstruction of the real object 
trajectory from the trend. As far as we know, this method has never been used before in the 
metocean field. 

In order to check the validity of the proposed methodology, we consider a situation with deter- 
ministic currents and negligible wave effects and we study a simplified model where the acceleration 
of the floating object is equal to the sum of the acceleration of the water (due to the wave tide and a 
perturbation of smaller order which represents other water currents) and the difference between the 
wind speed and the object velocity. This model is described more precisely in the second section. 
Then, in the third section, we present the basic ideas of the Averaging Method and we compare 
computational cost of this method with the one which consists in integrating the real trajectory 
directly. In Section 4, we validate our method using numerical experiments. First, the models of the 
meteorological forcing fields (tide wave, perturbation and wind) are specified. Then, the accuracy 
of the Averaging Method is checked on the basis of numerical comparisons. Finally, in Section 5 the 
methodology is illustrated with an example in which the probability of running aground is estimated 
for an academic domain. 

2 Model 

The model on which we shall implement the method evoked above is a very simplified model for large 
time drift in ocean, above continental shelf in strong tide zone, of an almost completely submerged 
object submitted to wind. This model is extracted from exhaustive scaling analysis for large time 
floating object drift which we shall present in a forthcoming paper. 

In short, the evolution of the position X(t) = X(f; x, v) G and the velocity V(t) = V(i; x, v) 
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e of the considered object, having x and v as initial position and velocity, is given by: 



dt " dt 



dm 



m(i, X)) + A(w(t, X) - V) = —{t, X) + (Vm(i, X)) V + A(w(t, X) - V), (2.2) 



where m = m.{t, x) is the ocean velocity field and w = w(t, x) the wind velocity field. Vm stands 
for the Jacobian matrix of m. This equation says nothing but that the object is submitted to the 
sea water acceleration and to the wind force quantified by a constant A. 

The time scale t on which wc want to observe the drift phenomenon is about 3 months and the 
object is submitted to tide oscillations whose period T is about 12.5 hours. Then a small parameter 
e appears in our problem: the ratio £ tide period on observation time scale. The magnitude of e is 
about 1/200. Wc consider that the velocity measurement of the object, wind and water arc all done 
with the same unit: v. The observation length scale I, which has to be the characteristic length 
of continental shelf, is about several hundred kilometers. This length has to be compared with the 
characteristic distance Tv the water covers during a tide period, which is of some kilometers. It 
seems then reasonable to consider that the ratio ^2 is also about e. 

Having those scale considerations at hand, we introduce the following rescaled variables t' , x' 
and v' expressing time, position and velocity in unit t, I and v 

t = tt', X = Ix' and v = vV, (2.3) 

the rescaled trajectory (X'(i'; x', v'), V'(i'; x', v')) defined by 

IX'(t';x',v') = X(tt';Ix', w'), vV'{t';x' = V(«';Ix',tJv'), (2.4) 

and the rescaled fields m' and w' defined by 

vin'{t',x') = m(tt',Ix'), vw'{t',x') = w(tt',Ix'). (2.5) 

We consider that the sea velocity writes 

m'(i', x') = M{t', L', x') + £N(i', L', x') = M{t', ^, x') + £N(i', ^, x'), (2.6) 

where M(i', j,x') G is the rescaled sea water velocity exclusively due to the tide wave. It 
is supposed to be regular. Considering its dependency with respect to the oscillating time vari- 
able, we suppose that 6' ^ M(f',0',x') is a 1— periodic function satisfying ^ M(f',6',x')d^' = 

where j M{t' ,9' ,x') d0' = J M{t' ,9' ,x') dO' . The field sN{t' ,'-,x'), where N{t',9',x') is also 

1— periodic in 6', is the sea water velocity perturbation induced by meteorological factors. As con- 
cerns the wind velocity field, we consider that w' also involves two time scales, i.e., 

w'(i',x')=W(i',^,x'). (2.7) 



Nonetheless, wind time series bring out as unrealistic considering a periodic dependency of W(t', 9' , x') 
with respect to 9' . In practice, wc only consider that, for any t' and x', it admits an average value 

W(t', 6' , x') d9' which actual definition is discussed later on. 

Finally, we may deduce the equation satisfied by the rescaled trajectory and involving the rescaled 
fields. Removing the ', this equation, which is the model on which we shall implement our method. 



/ 
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reads: 



(2.9) 

Notice that system p.8p - p.9p is rescaled, and that, in it, every variable and field characteristic scale 
is of order 1 . 

Clearly, the model under consideration is too simplistic to be used for operational applications. 
Nevertheless, it contains most of the physical ingredients of drift of object in the ocean: joint action 
of sea and wind, two time scales, possibility of using not so unrealistic sea velocity fields. Moreover, 
it seems to be relatively straightforward to incorporate a realistic sea velocity in it, with a tide period 
that weakly evolves with time and to apply it in a real geographical geometry. Hence, the validity 
of the methodology we present in this paper is not limited by the simplifications we consider. 



3 Asymptotic analysis 

Having the goal of using the Monte Carlo Method to estimate the probability of events linked to 
the trajectory of the drifting object in mind, we need to compute the object trajectory for a large 
number of wind conditions. The solution (X, V) contains i— frequency oscillations. Then, solving 
(|2.8p - (|2.9p directly by numerical methods forces the use of a very small time step. For instance, if 
the explicit Euler scheme is used, since the characteristic size of the left hand side of (|2.9|) and of 
its gradient are about - and the size of its time derivative is about , the classical error estimate 
yields, for a time step At small enough, an error about 

At/ At\^ At 1 , , 

— 1 + — ^ — e^. (3.1) 



e \ e _ 

Hence, if we want to obtain a precision about a time step At about e^e~i is needed. This is 
really too small for operational applications. Hence we shall write an expansion of (X, V) and find 
non oscillating equations satisfied by the terms of this expansion. This way, the previously evoked 
constraint imposed on the time step vanishes. 

It is an easy game to see that system (|2.8p - (|2.9p enters the framework of an oscillatory-singularly 
perturbed dynamical system 

4(^1 =a(t,-,X,V)+eai(t,-,X,V) + ib(t,-,X,V), (3.2) 
at \y j £ e e e 

very close to the one studied by Frenod [7] which originates from gyrokinetic plasma questions (see 
also Poincare [T^ , Krylov and Bogoliubov [TJ , Bogoliubov and Mitropolsky [3] , Sanders and Verhulst 
[TO] . Schochet [201, Frenod and Sonnendriicker [SKTUIIII], Frenod and Watbled [Hj and Frenod, 
Raviart and Sonnendriicker [S] for presentations of methods allowing to remove time oscillations). 
With very little changes we may apply the results it contains to deduce that 

X(i)=XO(t,-)+£Xi(<, -) + ..., V{t)=V"{t,i)+eV\t,i) + ..., (3.3) 

£ e e £ 
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where oscillating functions X°, V°, and are linked to non oscillating functions Y", U°, 
and by 

XO(t,^)=YO(i), (3.4) 
Y°{t, 6) = M(t, e, Y\t)) + U°(t), (3.5) 

and 

X^{t,e)=Y^{t)+ M{t,a,Y°{t))da, (3.6) 
Jo 

V\t,e) = {WM{t,e,Y''m{Y\t)+ [ M(t,<7,Y0(t))d<7} 

+ (i) + N{t, e, Y°(i)) - N(t, 0, Y°{t)) + J (W{t, a, Y° {t)) - j W{t, Y°{t)) d^) da 
1-0 

- / M{t,a,Y^{t))da. 
Jo 

(3.7) 

Then Y°, U°, Y^ and are the solution to 

ciYO 



= U°, (3.8) 

=J^wit,e,Y°)de-v°, (3.9) 



dt 
dlJ° 



dY 

~di 



-= j {VM(i, 6, Y°)}{ M(t, (7, YO) rfo-} d6l 

+ + ^ N(t, 0, Y°) d0 - N(i, 0, Y°) + j (W(<, a, Y^) - ^ W(<, ^, Y°) d^) da d0 

-jj M{t,a,Y°)dade-[J^ J VM(t, a, y^) darf^}{UO} - ^ ^ ^(t,ay)dade, 

(3.10) 



- = { _^ VW(t, ^, Y°)de}{Yi} + ^ {VW{t, e,Y°)}{J^ M{t, a, Y°) dajdO 

- { ^ {VM{t, e, Y°)}{ M{t, a, Y°) dc7}d6' + + ^ N(t, 6*, Y°) ^6* - N{t, 0, Y°) 

+ 'W{t,a,Y°)da-eJ'W{t,a,Y°)da)d9-J' m{t,a,Y°) dade"^ 

+ {vN(i,0,YO) 

~f^^lo ^(■'^'■)^^~^/ ^(■''^'■)'^°')(*'^°)'^^ + / ^ M(i,c7,Y°)rfad6l}{U°} 

d( ['w{;a,-)da-eI'W{;a,-)da) 
+ ^(M,Y°)-^ ^^^2 _J i(i,YO)d^ 

+ f J -^{t,a,Y°)dade, 
(3.11) 



equipped with initial conditions Y°(0) = x, U°(0) = v, Yi(0) = and \J\0) = 0. 

As said previously, this result may be inferred from [J. It may also be deduced by expanding 
fields and functions in (j2.8p and (|2.9p in a correct way. This is done formally in the Appendix. 

Equations l|3.4p and (|3.5p mean that the order trajectory does not oscillate and that the order 
velocity is the tide velocity added with a non oscillating velocity U" which is generated by wind. 
Wind acts on U° only through its averaged value. This is translated in (|3.9p . By the way, notice that 
since the average value of M is 0, equations (|3.8p - (|3.9p only involve the average wind. As concerns 
order 1 terms, the situation is more complex. First since M(<, a, Y'^(t)) da may be interpreted as 
the position of a sea water particle placed in Y"(i) at the beginning of a tide cycle (9 = 0), equation 
(|3.6p means that the order 1 position is this water particle position plus a non oscillating function 
Y^{t). Regarding the terms (|3.7p contains, the first one describes the way the space variation of the 
tide velocity acts, the second one is the non oscillating part of the velocity. The third and fourth terms 
quantify the action of the sea velocity perturbation. Concerning the next term, we need to remember 
that the action of the averaged value of the wind velocity is taken into account in order equation 
()3.9p . Then we notice that W(i, cr, Y'^(i)) ~ / W(t, Y°(f)) quantifies the wind action around its 

averaged value at each time of a tide cycle. Hence (^W(i, cr, Y°(t)) — / W(t, Y"(t))) d<j^ da is 

the cumulated action of the wind around its averaged value. This quantity acts at order 1. The 
last term of p.7p can be interpreted as the previous one recalling that the mean value of the sea 
velocity is 0. It is hard to give intuitive explanations for the evolution equations (|3.10p - (|3.1ip . We 
only notice that they involve mean value of non-linear interactions between fields which quantify the 
mean joint action of sea and wind. This non intuitive quantification is made possible thanks to the 
asymptotic analysis presented in 7 and in the Appendix. 

The characteristic size of the right hand sides in p.Sp - p.lip are 1. Hence, to get the same 
precision about e^, with the same Eulcr scheme using p.3p - p.lip to compute (X,V), in place of 
using (I2.8p - ()2.9p as evoked in the beginning of this section, the needed time step At is about 
Comparing this with the value e^e~^ found in the beginning of the section, despite the heavy form 
of ()3.8p - (l3.1ip . this is an appreciable gain even more so given that our probabilistic forecast method 
needs numbers of simulations for numbers of wind time series. 



Concerning the validity of the expansion (|3.3p . if W(t, 6, x) was a regular function, periodic with 
respect to 9, asymptotic expansion p.3p could be rigorously justified, applying Jj, by the following 
inequality 

sup ||X(i)-XO(t,-)|| <ce, sup ||X(t)-X"(t,-)-eXi(t,-)|| <c£2, (3.12) 
te[o,i] ^ te[o,i] s ^ 

sup I|V(<)-V°(i,-)|| <ce, sup ||V(t)-VO(i,-)-eVi(i,-)|| <c£^ (3.13) 
te[o,i] ^ t6[o,i] ^ ^ 

which would be true for e small enough and for a constant c independent of e, where || • || stands for 
the Euclidean norm in M.^. 



As we cannot consider the wind time series to have this form, we will estimate the error using 
numerical experiments. This is done in the next section: several realistic wind time series are 
simulated using a method described in Monbet et al. [13] and the corresponding error is calculated 
for each of them. For this, we have to decide how the average values are to be computed. If the way 
is clear for fields linked to M or N, for which we simply take ^ d9 = d9, the situation is more 
obscure for wind linked fields. For them, we will take as averaged value at time t the mean value on 
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the interval centered in t and with length p, i.e. 



f 1 s 

/ W(t,6,x.)dd ^ - W(s, -,x)ds, (3.14) 

J P Jt-p/2 £ 



for a parameter p that has to be adjusted experimentally. As concerns the integrals in 9 of the wind 
linked fields, they have to be replaced by integrals respecting Vi(t, |) = Ui(t) when | is an integer, 
which is a constraint imposed by the Averaging Method. Hence, we replace 

/ W(i,(T,x)dcr) by / W(s, -,x)ds, (3.15) 

Jo J Je[^] e 

where [-] stands for the integer part of -. Those choices have also to be experimentally validated. 



4 Numerical validation 

We begin this section by introducing the sea velocity fields and briefly presenting the stochastic 
method used to generate the wind time series. Then, the expansion p.3p is validated by numerical 
experiments. More precisely, several wind time series are generated. Then, for each of them, 
we compute the real trajectory using (|2.8p - (|2.9p and compare it to the trajectory obtained with 
expansion (|3.3p and the following equations. 

4.1 Metocean fields 

Let us first describe the metocean fields that have been used in our numerical experiments. Regarding 
the rescaled ocean velocity induced by the tide wave and its perturbation, we have chosen the 
following simple parametric fields: 

/ sin(27r6') + isin(47r6, , 
M(t,0,x) = (2 + sin(6^t)) .Ti ^ ^ 4 V ,^ ^^^^^ 

y isin(27re' 
/ sin(27r6') ' 

N(i,6',x) = (2 + cos(67rt)) X2 , (4.2) 

ysin(27r6') ^ 

where xi and X2 are the first and second components of x. 

/ sin(27r6') + isin(47r6^ 

The trajectories associated with the field ' I are non circular loops which 

y i sin(27r6') 

remind those we can observe in Nihoul [16j or Salomon and Breton [18]. The velocity field M is 
this simple vector field modulated by a time and position dependency in order to see the influence 
of time derivative and gradient on the object trajectory. N is also a simple field which gradient is 
orthogonal to the one of M. 

In order to simulate realistic wind time series, we have assumed that they can be decomposed as 
the sum of two components, e.g W(t, |, x) = Witit, x) + Wst(|, x) where 

• Wit(<,x) represents the wind evolution at a synoptic-scale, e.g at the scale of the high- and 
low-pressure systems. The typical dimension of these systems ranges approximately from 
1000km and 2500km and their duration is a couple of days to at most a couple of weeks. 



Wsi(|, x) represents the small scale evolution of the wind (e.g mesoscale and microscale winds). 
This scale includes phenomena such as thunderstorms, squall lines, land and sea breezes... 
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Such a decomposition is discussed more precisely by Breckling (see Breckling |3] and [T]). Dif- 
ferent methods have been proposed in the Hterature to simulate realistic wind time-series at the 
synoptic scale (see Monbet, Ailliot and Prevosto [TS] and references therein). In this study, we have 
first assumed that the wind Wit(i, x) is homogeneous in space; i.e. Wnit, x) = Wi((i) for all x and 
t, what seems realistic according the size of the domain which is supposed to be about some hundred 
kilometers. Then, we have used a non-parametric resampling method to simulate the process 'Wnit), 
namely the Local Grid Bootstrap algorithm proposed by Monbet, Ailliot and Prevosto [15]. This 
method has already been validated on various datasets, and it was found that it can successfully be 
used to simulate realistic wind time series. In the present study, it has been calibrated on a dataset 
which describes the wind condition during the summer at a point of coordinates (46.25 A^, 1.67 W), 
located near the French Atlantic coast. It describes the synoptic wind conditions during the last 20 
years recorded every Ati = 6 hours. 

Then, in order to simulate the small scale variations, we have used an Autoregressive model. 
More precisely, for simplicity reasons, we have first assumed that this field is homogeneous in space. 
This assumption is unrealistic and could be refined later on. Then, we have assumed that, for 
k e N*, 

WstikAh) = aWstiik - l)At2) + E{kAt2) (4.3) 

where {E{kAt2)}k£K-' denotes a zero-mean Gaussian white noise with covariance matrix a^l2- In 
practice, we have simulated this small-scale component with a time-step At2 = e/100. This consists 
in giving the wind value every 4 min approximately, and the parameters a and a have been chosen 
such that the process W^t has a memory of a few hours and such that the standard deviation of its 
marginal distribution represents approximately 10% of the one of the process W. 

Finally, the procedure described above makes it possible to simulate the processes Wit(t) and 
Wji(-) for t e fcAti and - G kAt2, respectively. The values of these processes for other values of t 
have been calculated by linear extrapolation. An example of simulated wind time-series is given in 
Figure [H 

4.2 Numerical results 

In order to validate the asymptotic expansion given in the previous section, we have computed the 
solutions of (I2.8|) - (|2.9|) and of ()3.3p - (l3.1H) for N — 100 artificial wind time series. Let us denote 
Xi, Vi, X°, V°, Xj, "Vj, for i G {1, ...,N}, the corresponding numerical approximations of X, V, 
X", V°, X^ and V^, respectively. We have tested different algorithms to solve these DDEs, and 
the best results have been obtained with an explicit Runge-Kutta (4,5) formula. In practice, we 
have used the Maltlab's function ode45 which is based on an explicit Runge-Kutta (4,5) formula, 
the Dormand-Prince pair [6|. In the numerical results given hereafter, we have used e — 1/50. This 
choice makes it possible to compute the solutions of the system (j2.8|) - (|2.9|) with a good precision in 
a reasonable computational time, and permits also easier graphics representations. We have used 
Xi — (1, 1) and Vi = (0,0) as initial values. 

In Table [Tl the norms of error in object position and velocity for order and 1 expansions are 
given. Let us discuss more precisely the results obtained with e = 1/50. It shows that the asymptotic 
expansions obtained with p = e/10 (corresponding to an interval of 40 min), p = s/2 (corresponding 
to an interval of 6hl5min) and p = e (corresponding to an interval of 12h30min) are about 5£^ 
worth, and close to each other. For comparison, we also compute the solution to system (I3.3p - p.lip 
when the wind is null, and we also found an error equal to 0.0022 « 5e^. For p = Ae (corresponding 
to an interval of 2 days), the error is significantly higher. As concerns the computational coast, 
it decreases as p increases, so that a good compromise seems to use e/2 < p < e; i.e. the wind 
at a synoptic scale. Such a value of p permits to compute the solutions of the system p.3|) - (|3.1H) 
with a good precision in a computational time significantly lower than the one corresponding to the 
system (|2.8p - (|2.9p . For instance, if an Eulcr scheme is used, the exact system (|2.8p - (|2.9p requires 
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p = 


6 

10 


0.1115 [0.1115,0.1116] 


0.0460 [0.0454,0.0467] 


0.0174 [0.0173,0.0174] 


0.0023 [0.0021,0.0029] 


p = 


2 


0.1118 [0.1113,0.1134] 


0.0460 [0.0448,0.0473] 


0.0174 [0.0173,0.0174] 


0.0024 [0.0021,0.0032] 


p = 


e 


0.1125 [0.1108,0.1168] 


0.0460 [0.0441,0.0488] 


0.0173 [0.0171,0.0174] 


0.0026 [0.0020,0.0047] 


p = 


4e 


0.1192 [0.1102,0.1280] 


0.0542 [0.0438,0.0691] 


0.0174 [0.0166,0.0196] 


0.0065 [0.0024,0.0144] 


W E 


E 


0.1115 


0.0096 


0.0174 


0.0022 



Table 1: Mean value, minimum and maximum values (mean [min, max]) of the errors 
suptg[o,i] ||V,(t) - V^{t)\\ (second column), suptg[o,i] l|V»(t) - V,"(t) - eV}it)\\ (third column), 
suptg[o_i] ||X,(t) - X.°{t)\\ (forth column) and sup^gfo.!] W^ti*) - - £^}{t)\\ (fifth colmim) for 

different values of p. The last line for e = gives the error for a zero wind field. 



about 1000 more iterations than the approximate system (j3.3|) - p.lip to achieve the same accuracy. 
And, for the problem considered in this paper, an iteration of the exact system is 20 less expensive 
in terms of computational time. But this last remark is not general since the computational time 
highly depends on the nature of the tide and current fields M and N : here they are modeled by a 
quite simple analytical formula. 

In Figure [Tl we have plotted the object trajectory associated with the wind time series shown on 
the top of Figure m using p = e/2. More precisely, the sohd hne represent X computed by directly 
solving (j2.8p - (|2.9p . The dashed line represents the average trajectory Y*^ + eY-'^ obtained solving 
p.8p and (|3.10p . We can see that this averaged trajectory follows nicely the trend of X. Then 
reconstructing X° and X^ using (|3.4p and (|3.6p . we have represented the trajectory X" + eX^ in 
dotted line. The superimposition with X seems to be almost perfect. In order to go further in the 
result analysis, we turn to Figure [H On the top of this figure, we have shown the wind as a function 
of time. The second figure represents, in dashed line, the average trajectory of the position of the 
first component and, in solid line, the first component of the trajectory itself. We can see on this 
trajectory not only the long term trend but also the two time periodic phenomena in game: tide 
oscillations (rapid oscillations) and tide coefficient amplitude (modulated amplitude). Finally, the 
last plot exhibit e^— order error (here e = 1/50 = 2. 10~^, then 10~^ = 2.5£^). Moreover, we can 
see that this error function is a periodic function with modulated amplitude. This spurs us thinking 
that we could improve the accuracy of the reconstructed trajectory, if it is needed, considering the 
next terms in the expansion of X and V. Figure [3] shows the first component of the average wind, 
the average velocity (dashed line) and the velocity itself (solid line). The third plot shows the error 
on the velocity. This figure permits to visualize the action of the wind on the averaged velocity, 
which first increases and then decreases. 

In order to make our numerical validation more convincing, we also made simulations for varying 
e. In Table m the norms of the error on object position and velocity for order and 1 expansions 
are given for e = 1/100 and e = 1/25 . The errors are proportional to e for order and to for 
order 1 as it was expected from the theory. 
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Figure 2: Example of object's trajectory (time series plot) (Top: zonal wind first component, Second: 
first component of zonal object position (Solid line: X, Dashed line: Y° +eY^), Third: zonal error 
||X-X" -eXHl). £ = l/50,p = e/2 



11 



3 - 




Figure 3: Example of object's speed (time series plot) (Top: smoothed zonal wind component (p — 
e/2), Second: first component of zonal object's speed (Solid line: V, Dashed line: U°+eU^), Third: 
zonal error || V - - eV^ ||). e = l/50,p = e/2 
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1 

100 


Speed (order 0) 


Speed (order 1) 


Position (order 0) 


Position (order 1) 


■ = 


e 

10 


U.UOO/ [U.UOO(,U.UOO(J 


n 0909 rn 0970 n 0*^^81 


n nnsfi rn nns"^ n nnsf^i 

U.UUoU [u.UUoO,U.UUOOJ 


n nni i r n nni i n nni a\ 


P = 


e 
2 


U.UOOo [U.UOO^,U.UODoJ 


n n9QA rn 0970 n c\'\f\0\ 


n nnRfi rn nns"^ n nnsf^i 


n nni 9 rn nni n n nni ^1 


P = 


e 


0.0560 [0.0552,0.0574] 


0.0296 [0.0273,0.0363] 


0.0085 [0.0084,0.0086] 


0.0013 [0.0010,0.0020] 


P = 


4e 


0.0595 [0.0552,0.680] 


0.0355 [0.0304,0.0516] 


0.0086 [0.0082,0.0098] 


0.0026 [0.0010,0.0072] 


e = 


1 

25 


Speed (order 0) 


Speed (order 1) 


Position (order 0) 


Position (order 1) 


P = 


e 

10 


0.2254 [0.2249, 0.2257] 


0.1143 [0.1022 ,0.1418] 


0.0358 [0.0358,0.0359] 


0.0046 [0.0040,0.0057] 


P = 


2 


0.2254 [0.2226,0.2277] 


0.1142 [0.1017,0.1410] 


0.0358 [0.0357,0.0359] 


0.0045 [0.0039,0.0058] 


P = 




0.2257 [0.2201,0.2310] 


0.1142 [0.1011,0.1379] 


0.0358 [0.0356,0.0360] 


0.0044 [0.0037,0.0069] 



Table 2: Mean value, minimum and maximum values (mean [min, max]) of the errors 
suptg[o,i] ||V,(t) - W°{t)\\ (second column), sup^gfo ij ||V,(t) - VO(t) - eVKOII (third column), 

suptg[o,i] l|X»(i) - X^(i)|| (forth column) and suptg[o,i] \\%{t) - X^(t) - e±l{t)\\ (fifth column) for 
different values of p. 



5 Application: long term forecast of an object's drift with 
Monte Carlo Method 

In the previous section, we have shown that the expansion of the trajectories of an object in near 
coastal ocean submitted to wind makes it possible for us to compute quickly good approximation 
of this trajectory. It is then possible to make this computation for a large number of synthetic 
wind time series, and thus deduce the probability of a given scenario (this may be the probability 
of presence in a given area, of collision or of running aground, for example). 

As an example, in this section we compute a running aground probability. More precisely, we 
have fixed an arbitrary coast line, e.g. the circle of center (1, 1) and radius 0.3 and we have focused 
on the running aground of the object on this coast. We have generated one thousand wind time 
series with the same stochastic model as in the previous section and have also used the same sea 
velocity fields. We have computed the thousand associated trajectory. For this, we have used the 
method described in the previous sections, based on asymptotic expansion (j3.3p . The methodology 
is illustrated in Figure [3] where one hundred of the thousand trajectories are drawn. Proceeding this 
way, we have obtained enough trajectories to get reliable estimates of the quantities of interest. First 
we can deduce an estimation of the probability of running aground by calculating the percentage of 
trajectories that reach the circle in the time interval [0, 1]. We found 64.3%. Then we can estimate 
the probability that the object runs aground in a given area. The obtained results are shown in 
Figure O On the left, we have plotted a wind rose which shows the joint repartition of direction 
and intensity. The length of each bar is proportional to the proportion of wind having a direction 
around the direction the bar indicates. The bars are then shared into four subbars which lengths are 
proportional to the proportion of wind having the corresponding magnitude. We can see that the 
wind is generally blowing from the north. The second histogram must be read in the following way. 
The length of each bar is proportional to the proportion of trajectories that run aground around the 
point of the coast which is in the bar direction. We then see that the proportion of running aground 
trajectories has a maximum around the South-South- West of the coast line. 
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Figure 4: Object's trajectory corresponding to 100 synthetic wind time series and corresponding 
running aground points(*) 
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6 Conclusion 



In this paper, we have explored a method to make probabilistic forecasts of long term drift, in 
near coastal ocean, of object submitted to tide and wind. The method is based on generating a 
large number of realistic wind time series and, for each of them, to compute the associated object 
trajectory. The probabilities of trajectory linked events may then be estimated using this large 
number of trajectory realizations. In order to achieve this goal, we have used a method proposed by 
Monbet, Ailliot and Prevosto [15] to generate realistic wind time series and the Averaging Method 
based on an asymptotic expansion of trajectory of Frenod [7] to remove tide oscillations and compute 
quickly the object's trajectory. This method has allowed us to compute the probability for an object 
to run aground on an academic domain. 

The quality of the results are good enough for us to contemplate going further in exploring 
this new methodology. Among the things to do in order to achieve the objective of computing 
operationally events probabilities in real coastal areas we wish to broach the following. 

First, we will present a complete scale analysis of a coupled system Shallow Water Equations - 
Newton Principle modeling the joint running of ocean and drifting object. The aim is to identify 
several regimes under interest such as "storm regime in coastal zone" or "stillness above continental 
shelf" and so on, and the corresponding equations. This would give a way to compute the sea 
velocity fields M and N, M being exclusively due to tide wave and N to meteorological factors. 
Notice also that N could be divided into two parts: one linked to pressure variations and another 
linked to wind. Those aspects would then have to be incorporated in a software in order to access 
the fields M and N. Then we shall study the abilities of the method presented in the present paper 
to fit to pseudo-periodic fields M and N. We are optimistic because the method already works 
well with pseudo-periodic wind fields. The realism of the synthetic wind time series used as input 
of the numerical model could also be improved. In particular, the space-time model proposed by 
Ailliot, Monbet and Prevosto [2] could be used to simulate non-homogeneous wind fields. Finally, 
the method could be implemented on real coastal areas. 
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A Appendix : Computations leading ( 13.41) - ( 13.111) 



In this Appendix we formally lead the asymptotic expansion that gives equations p.4p - (|3.1ip . 
First, if we define Y(t) and U(i) so that 

X{t)=Y{t), (A.l) 

V(<) = M(i,J,Y(i)) + U(i), (A.2) 

and inserting those in (|2.8p and (|2.9p . we deduce 
dY f 

— =M(t,-,Y)+U, (A.3) 

dt £ 

+W(i,-,Y)-M(t,i,Y)-U. 

£ e 
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Then, we assume that Y and U may be expanded in the foUowing way: 

Y(t) - Y"(i) + e(Yi(t) + A\t, -,Y"m + e{Y^{t) + A\t, ^ ,Y\t) ,Y\t))) + • • • , (A.5) 

£ e 

U(i) = UO(i) + £(Ui(i) + Bi(i, -, YO(t))) + e{\]^{t) + B2(t, YO(t), Yi(i))) + • ■ • , (A.6) 

£ £ 

for functions A''(i, 6*, Y^), (i, 6*, Y^, Y^), B°{t,6,Y°) and Bi(i, 6i, Y", Y^) being 1- periodic with 
respect to to be defined latter. Using those asymptotic expansions in (lA.ip and (jA.2p yields (13. 4|) 
and (|3.5p . Using again (jA.Sp and (|A.6p in (jA.Sp and (|A.4[) . expanding every functions using a Taylor 
expansion yields 

^ + =^ + =^<^") + ^<^°' + ^('^' <^°') (^) + ^1^<^°' ^' ' + ^ ■ ■ • 

= M(Y*') + e(VM(Y°)) (Y^ + A1(Y")) + U" + £{\J^ + 1A\Y°)) + ■■■ (A.7) 



= |^(Y") + e(f:|^(YO)) (Yi + Ai(YO)) + £^(Y0) + e(VN(YO)) (M(Y0) + U") 

+ W(Y") + e(VW(Y°)) (Y^ + A\Y^)) - M(Y") - e(VM(Y°)) (Y^ + A\Y'^)) 

-U"-£(Ui+Bi(Y")) + --- , (A.8) 

where + • • ■ contains every terms being of order greater than 1 in £. 
Identifying the terms of order in £ gives 

^ + ^(YO) = M(Y°) + UO, (A.9) 

^ + = ^(Y«) + W(Y°) - M(YO) - U". (A.IO) 

Integrating those equations with respect to 9 from to 1, we obtain (|3.8p and (|3.9p . Once this is 
done, equations (|A.9p and (jA.lOp yield 

A^{t,0,Y")= M{t,a,Y°{t))dcr, (A.ll) 
Jo 

B^t, e, Y") - V\t) + N(t, 0, Y"(t)) - N(t, 0, Y°(t)) 

+ y (wit,a,Y°it))- Wit,i,Y°{t))d^jda~ J M{t, a,Y° (t)) da, 

and then dSS]) and ((37)) . 

Looking now at the terms of order 1 in £ in (jA.7p and (|A.8p . after replacing A-^, B^, dY°/dt and 
dXJ'^/dt by their expressions, results in two equations. Those equations, when integrated in 0, lead 
in a heavy but straightforward way (|3.10p and (|3.1ip . 
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